rm(list=ls())

root = "/Users/Tom/Documents/YNP_Bison_Model/A_final_model_runs_for_revised_manuscript/"
root2="/Users/Tom/Documents/YNP_Bison_Manuscript/A_Ecological_Monograph_Submission/Revised manuscript/Results/"


path = function (stem,r=root){
return(paste(r,stem, sep=""))	
}

library(rjags)
library(coda)
library(xtable)

setwd(root)

#values in text

#Get median of z

load(path("Freq and Dens Dependent/One_beta_jags.Rdata"))
 

 z.med=summary(One_beta.jags$z,quantile, .5)$stat

 save(z.med, file=path("median of z.Rdata", r=root2))


#TABLES==================================
#Table for model comparisons=============
#load(path("Model and analysis code/Frequency Dependent/One_beta_jags.Rdata"))
load(path("Frequency Dependent/One_beta_jags.Rdata"))
##"))
#####DIC for model comparisons
freq.deviance = summary(One_beta.jags$deviance,mean)$stat
freq.deviance.var=mean(summary(One_beta.jags$deviance,var)$stat)
freq.dic=freq.deviance+1/2*(freq.deviance.var)
freq.mspe=summary(One_beta.jags$mspe2,mean)$stat
freq.P_yphi_p = summary(One_beta.jags$P_yphi_p,mean)$stat

load(path("Density Dependent/One_beta_jags.Rdata"))
dens.deviance = summary(One_beta.jags$deviance,mean)$stat
dens.deviance.var=mean(summary(One_beta.jags$deviance,var)$stat)
dens.dic=dens.deviance+1/2*(dens.deviance.var)
dens.mspe=summary(One_beta.jags$mspe2,mean)$stat
dens.P_yphi_p = summary(One_beta.jags$P_yphi_p,mean)$stat

load(path("Freq and Dens Dependent/One_beta_jags.Rdata"))
comb.deviance = summary(One_beta.jags$deviance,mean)$stat
comb.deviance.var=mean(summary(One_beta.jags$deviance,var)$stat)
comb.dic=comb.deviance+1/2*(comb.deviance.var)
comb.mspe=summary(One_beta.jags$mspe2,mean)$stat
comb.P_yphi_p = summary(One_beta.jags$P_yphi_p,mean)$stat



DIC.table=matrix(nrow=3,ncol=4)
rownames(DIC.table)=c("Frequency","Combined", "Density")
colnames(DIC.table) = c("Deviance", "DIC", "MSPE $\\phi$", "LPPD$\\phi$")
DIC.table[2,]=c(comb.deviance, comb.dic, comb.mspe, log(comb.P_yphi_p))
DIC.table[1,]=c(freq.deviance, freq.dic, freq.mspe, log(freq.P_yphi_p))
DIC.table[3,]=c(dens.deviance, dens.dic, dens.mspe, log(dens.P_yphi_p))

print(xtable(DIC.table[,c(1,3,4)],align=c("l","c","c","c"),digits=3, display=c("s","fg","fg","fg")), NA.string="NA", floating=FALSE,sanitize.colnames.function = function(x) {x})
save(DIC.table, file=path("DIC_table",r=root2))


#Evaluation of convergence
load(path("Frequency Dependent/One_beta_coda.Rdata"))
gelman.diag(One_beta.coda)
heidel.diag(One_beta.coda)
acfplot(One_beta.coda, lag=100)

#Evaluation of convergence
load(path("Density Dependent/One_beta_coda.Rdata"))
gelman.diag(One_beta.coda)
heidel.diag(One_beta.coda)
acfplot(One_beta.coda, lag=100)


#Evaluation of convergence
load(path("Freq and Dens Dependent/One_beta_coda.Rdata"))
gelman.diag(One_beta.coda)
heidel.diag(One_beta.coda)
acfplot(One_beta.coda, lag=100)

#Parameter estimates============================
load(path("Frequency Dependent/One_beta_coda.Rdata"))
load(path("Frequency Dependent/One_beta_jags.Rdata"))
zj=One_beta.jags
zc=One_beta.coda
#params=c("beta","b","m","psi","p","sigma.p","v")
#v=params
#jags.table=function(v,zj){
#	table = NULL
#	for(j in 1:length(v)){
#		a=summary(zj[[v[j]]],mean)$stat
#		b=summary(zj[[v[j]]],quantile,.5)$stat
#		c=colMeans(summary(zj[[v[j]]],sd)$stat)
#		d=summary(zj[[v[j]]],emp.hpd)$stat
#		add=cbind(a,b,c,t(d))
#		table=rbind(table,add)
#		}
#	colnames(table)=c("Mean", "Median", "SD", "HPDI", " ")
#	return(table)
#	
#}
#table = jags.table(v=params,zj=zj)
summary(zj$p,quantile,c(.025,.0975))$stat
table1=summary(zc)$stat[c(5,2:4, 7:12,14),1]
table2=summary(zc)$quantile[c(5,2:4, 7:12,14),3]


table3=summary(zc)$stat[c(5,2:4, 7:12,14),2]
table4=summary(zc)$quantile[c(5,2:4, 7:12,14),c(1,5)]
table=cbind(table1,table2,table3,table4)

colnames(table)=c("Mean", "Median", "SD","2.5%", "97.5%")
#The order is critical and it is critical that align and digits are within the xtable() functi, on.

rownames(table)[1:11]=c("$\\beta$","$f_{n}$","$f_{p}$","$f_{c}$", "$m$","$p_{1}$", "$p_{2}$", "$p_{3}$","$\\psi$", "$\\sigma_{p}$", "$v$")

print(xtable(table,align=c("l","r","r","r","r","r"),digits=2,display=c("s","fg","fg", "fg", "fg","fg")),NA.string="NA", floating=FALSE,sanitize.rownames.function = function(x) {x})
save(table, file=path("parameter_table.Rdata",r=root2))

#statistics for text
q=summary(One_beta.jags$diff.pos,quantile, c(.025,.5,.975))$stat
x=summary(One_beta.jags$diff.pos,mean)$stat
diff_pos=list(q=q,x=x)
save(diff_pos, file=path("diff.pos.Rdata",r=root2))

q=summary(One_beta.jags$diff.con,quantile, c(.025,.5,.975))$stat
x=summary(One_beta.jags$diff.con,mean)$stat
diff_con=list(q=q,x=x)
save(diff_con, file=path("diff.con.Rdata",r=root2))


#Posterior predictive checks================
ppc.table=as.data.frame(matrix(ncol=5,nrow=15))
ppc.table[,1]=c("Frequency dependent","","","","","Density Dependent","","","","","Combined","","","","")


ppc.table[,2]=c("Total population size", "Proportion juvenile","Juvenile serology", "Yearling serology", "Adult serology")


colnames(ppc.table)=c("Model","Data set","Discrepancy", "Mean", "SD")

load(path("Model and analysis code/Frequency Dependent/One_beta_JAGS.Rdata"))

m=matrix(nrow=3,ncol=3)
z=One_beta.jags
library(xtable)
ppc.table[1,3:5]=c(
summary(z$P.N.sq,mean)$stat,
summary(z$P.N.mean,mean)$stat,
summary(z$P.N.sd,mean)$stat
)
ppc.table[2,3:5]=c(
summary(z$P.ratio.calf.sq,mean)$stat,
summary(z$P.ratio.calf.mean,mean)$stat,
summary(z$P.ratio.calf.sd,mean)$stat
)
ppc.table[3:5,3:5]=c(
summary(z$P.sero.sq,mean)$stat,
summary(z$P.sero.mean,mean)$stat,
summary(z$P.sero.sd,mean)$stat
)
load(path("Density Dependent/One_beta_JAGS.Rdata"))
z=One_beta.jags
library(xtable)
ppc.table[6,3:5]=c(
summary(z$P.N.sq,mean)$stat,
summary(z$P.N.mean,mean)$stat,
summary(z$P.N.sd,mean)$stat
)
ppc.table[7,3:5]=c(
summary(z$P.ratio.calf.sq,mean)$stat,
summary(z$P.ratio.calf.mean,mean)$stat,
summary(z$P.ratio.calf.sd,mean)$stat
)
ppc.table[8:10,3:5]=c(
summary(z$P.sero.sq,mean)$stat,
summary(z$P.sero.mean,mean)$stat,
summary(z$P.sero.sd,mean)$stat
)
load(path("Freq and Dens Dependent/One_beta_JAGS.Rdata"))
z=One_beta.jags
library(xtable)
ppc.table[11,3:5]=c(
summary(z$P.N.sq,mean)$stat,
summary(z$P.N.mean,mean)$stat,
summary(z$P.N.sd,mean)$stat
)
ppc.table[12,3:5]=c(
summary(z$P.ratio.calf.sq,mean)$stat,
summary(z$P.ratio.calf.mean,mean)$stat,
summary(z$P.ratio.calf.sd,mean)$stat
)
ppc.table[13:15,3:5]=c(
summary(z$P.sero.sq,mean)$stat,
summary(z$P.sero.mean,mean)$stat,
summary(z$P.sero.sd,mean)$stat
)





print(xtable(ppc.table,align=c("l","l","l","c","c","c"),digits=2,display=c("s","s","s","fg","fg","fg")), NA.string="NA", floating=FALSE,sanitize.colnames.function = function(x) {x},include.rownames=FALSE)

save(ppc.table,file=path("Posterior pred check table.Rdata",r=root2) )



#Relative routes of infection

load(path("Frequency Dependent/One_beta_JAGS.Rdata"))
zj=One_beta.jags
p.infect = summary(zj$xNewIP,quantile, c(.025,.5,.975))$stat

p.table=matrix(0,nrow=8,ncol=3)
p.table[,1]=p.infect[2,c(1:5,7,8,6)]
p.table[,2]=p.infect[1,c(1:5,7,8,6)]
p.table[,3]=p.infect[3,c(1:5,7,8,6)]

rownames(p.table)= c("Vertical from sero-converter", "Vertical from recrudesced", "Horizontal juvenile", "Horizontal yearling", "Horizontal adult", "All vertical", "All horizontal", "Recrudesence")
colnames(p.table)=c("Median", "2.5%", "97.5%")

p.recrudesce= p.infect[,6]
save(p.recrudesce, file=path("Probability of recrudesence.Rdata", r=root2))

print(xtable(p.table,align=c("l","c","c","c"),digits=2,display=c("s","fg","fg","fg")), NA.string="NA", floating=FALSE, include.rownames=TRUE)
save(p.table,file=path("Routes of infection.Rdata",r=root2))

#Figures=====================
#functions, data, and constants needed for plotting below
load(path("Frequency Dependent/One_beta_jags.Rdata"))
load(path("2011 Bison Data with multinomial revised.Rdata"))
#general plotting function
plot.BCI=function(x,data.y, data.bar, model.y,y.high,y.low,xlab=years,ylab){
	y.errorbars(x,data.y,data.bar,y.low=y.low,y.high=y.high, 	xlab=xlab, ylab=ylab)

	lines (x,model.y[2,], typ='l', ylim=c(y.low,y.high)) #plots the median
	xx <-c(x, rev(x))
	yy <- c(c(model.y[3,]),rev(c(model.y[1,])))
	polygon(xx,yy,col=topo.colors(5,alpha=.3), border=NA)
}

y.errorbars = function(x,y,ybar, ylab=x, xlab=y, main=" ",y.low,y.high){
	plot(x,y,pch=16,,xlab=xlab, ylab = ylab, main=main, ylim=c(y.low,y.high), family="serif")
	arrows(x,y-ybar,x,y+ybar, code=3, angle=90, length=.01)
}

#Additional, redundant function, discovered after I wrote the one above!
plot_time_series = function(x, y, y.data, ybar, y.low,y.high, ylab, xlab, main){
y.errorbars(x=x,y=y.data,ybar=ybar,y.low=y.low,y.high=y.high, xlab=xlab, ylab=ylab, main="")
title(main=main, cex.main=1.25)
lines (x,y[2,], typ='l', ylim=c(y.low,y.high), xlab="Year", ylab=label) #plots the median
xx <-c(x,rev(x))
yy <- c(c(y[3,]),rev(c(y[1,])))
polygon(xx,yy,col=topo.colors(5,alpha=.3), border=NA)
}


#Set range for x=axis sequence
range=c(6:41)
#Set correpsonding range of years
yr=(1975:2010)

par(cex=1.25)

zm=One_beta.jags

#######end of data, constants, and functions needed for plotting


#plot posterior distributions=================
f=One_beta.jags
pdf(file=path("posterior density collage.pdf", r=root2), width =6, height = 8)
par(mfrow=c(4,3),cex.lab=1.5)
a=1

#beta
plot(density(f$beta, adjust=a), main = "A",  xlab = expression( beta), family="serif")
#plot(density(d$beta, adjust=a), main = "B", xlab = expression(paste("Transmission rate, ", beta," ", (year^{-1}))))
abline(h=dunif(5,0,10),lty="dashed")


plot(density(f$b[1,,], adjust=a), main="B", xlab=expression(italic(f[n])), ylim=c(0,12), family="serif", ylab="")
x=seq(0,1,.001); y=dbeta(x,77,18.1)
#lines(density(d$"b[1]", adjust=a), lty="dashed")
lines(x,y, lty="dashed")

#B[2]
plot(density(f$b[2,,], adjust=a),main="C", xlab=expression(italic(f[p])), ylab="", ylim=c(0,12), family="serif")
#lines(density(d$"b[2]", adjust=a), lty="dashed")
x=seq(0,1,.001); y=dbeta(x, 35.5, 20.5 )
lines(x,y, lty="dashed")
#B[3]
plot(density(f$b[3,,], adjust=a), main="D", xlab=expression(italic(f[c])), ylim=c(0,12), family="serif")
#lines(density(d$"b[3]", adjust=a), lty="dashed")
x=seq(0,1,.001); y=dbeta(x,3.2,11.3 )
lines(x,y, lty="dashed")

#m
plot(density(f$m, adjust=a), main = "E", xlab = expression(italic(m)), ylab="", ylim=c(0,12), family="serif")
#lines(density(d$"m", adjust=a), lty="dashed")
x=seq(0,1,.001);y=dbeta(x,49.5,49.5)
lines(x,y, lty="dashed")


#psi
plot(density(f$psi, adjust=a), main = "F", xlab = expression(psi), ylim=c(0,30), family="serif", ylab="")
x=seq(0,1,.001);y=dbeta(x,1.967,43.578)
lines(x,y,lty="dashed")

#p[1]
plot(density(f$p[1,,], adjust=a), main = "G", xlab = expression(italic(p[1])), ylim=c(0,20), family="serif")
abline(h=dunif(.5,.2,.95), lty="dashed")
#p[2]
plot(density(f$p[2,,], adjust=a), main = "H", xlab = expression(italic(p[2])), ylim=c(0,20),ylab="", family="serif")
x=seq(0,1,.001);y=dbeta(x,40,1.667)
lines(x,y, lty="dashed")
#p[3]
plot(density(f$p[3,,], adjust=a), main = "I", xlab = expression(italic(p[3])), ylim=c(0,20),ylab="", family="serif")
abline(h=dunif(.5,.2,.95), lty="dashed")

#sigma_p
plot(density(f$sigma.p[1,,], adjust=a), main = "J", xlab = expression(italic(sigma[p])), ylim=c(0,15), family="serif")
#lines(density(d$"p[2]", adjust=a), lty="dashed")
px=seq(.01,1,.01)
lines(px,dgamma(px,.001,.001), lty="dashed")


#v
plot(density(f$v, adjust=a), main = "K", xlab = expression(italic(v)), ylim=c(0,8), ylab="", family="serif")
#lines(density(d$v, adjust=a), lty="dashed")
x=seq(0,1,.001);y=dbeta(x,1.404,12.0578)
lines(x,y,lty="dashed")



dev.off()

#R_0
load(path("Analyses/R0.data.Rdata"))
 
pdf(file=path("R0.pdf",r=root2),height=6, width = 6)
par(mfrow=c(1,1))
hist(R0$f,breaks=100,freq=FALSE,xlab=expression(paste("Basic reproductive ratio, ", italic(R[0]))),xlim=c(1,3.5), ylim=c(0,2.5),main="",family="serif")
meanR0=mean(R0$f)
medianR0 = median(R0$f)
sdR0 = sd(R0$f)
ciR0=quantile(R0$f,c(.025,.975))
R0stats=list(meanR0,ciR0,medianR0)
save(R0stats, file=path("R0stats.Rdata"),r=root2)
text(3,2.0,paste("mean = ", signif(meanR0,3)), family="serif")
text(3,1.9,paste("median = ", signif(medianR0,3)), family="serif")
text(3,1.8,paste("SD = ", signif(sdR0,3)), family="serif")
text(3,1.7,paste("BCI = ", signif(ciR0[1],3),", ",signif(ciR0[2],3)), family="serif")


q_no_v=quantile(R0$no_v, c(.025,.5,.975))
q_no_psi = quantile(R0$no_psi,c(.025,.5,.975))
dev.off()


R0_list=list(meanR0=meanR0, medianR0=medianR0, ciR0=ciR0, q_no_v=q_no_v,q_no_psi=q_no_psi)
save(R0_list,file=path("R0.Rdata", r=root2))


R=data$removals$count
R90=R[R$year>=1990 & R$year <= 1999,]
bulls90=sum(R90$rm.bull)
cows90=sum(R90$rm.cow.adult + R90$rm.cow.yrlng)
bulls90/cows90

R00=R[R$year>=2000 & R$year <= 2010,]
bulls00=sum(R00$rm.bull)
cows00=sum(R00$rm.cow.adult+ R00$rm.cow.yrlng)
cows00/bulls00
zm=One_beta.jags
N1975=list(
q=summary(zm$N.total, quantile, c(.025,.5,.975))$stat[1:3,6],
mean=summary(zm$N.total,mean)$stat[6]
)
N2005=list(
q=summary(zm$N.total, quantile, c(.025,.5,.975))$stat[1:3,36],
mean=summary(zm$N.total,mean)$stat[36]
)
N2010=list(
q=summary(zm$N.total, quantile, c(.025,.5,.975))$stat[1:3,41],
mean=summary(zm$N.total,mean)$stat[41]
)
Ntotals=list(N1975=N1975,N2005=N2005, N2010=N2010)
save(Ntotals,file=path("Ntotals.Rdata",r=root2))




#lambda and demographics
load(path("Analyses/lambda_healthy.Rdata",r=root))
load(path("Analyses/lambda_infected.Rdata",r=root))

pdf(file=path("Lambda densities.pdf",r=root2),height=8, width = 7)
par(mfrow=c(2,1))
LH=lambda_healthy
hist(LH$lambda, breaks=100, main="A", xlab="", freq=FALSE,xlim=c(.90,1.25), family="serif")
lambda.mean=mean(LH$lambda)
lambda.median=median(LH$lambda)
lambda.sd=sd(LH$lambda)
lambda.ci=quantile(LH$lambda,c(.025,.975))
text(.96,15,paste("mean = ",signif(lambda.mean,3)),family="serif")
text(.96,12,paste("median = ",signif(lambda.median,3)),family="serif")
text(.96,9,paste("SD = ",signif(lambda.sd,3)),family="serif")
text(.96,6,paste("BCI = ",signif(lambda.ci[1],3), ", ",signif(lambda.ci[2],3)),family="serif")

adult.fem.median=quantile(LH$adult.fem,c(.5))
adult.fem.ci=quantile(LH$adult.fem,c(.025,.975))
yearling.fem.median=quantile(LH$yearling.fem,c(.5))
yearling.fem.ci=quantile(LH$yearling.fem,c(.025,.975))
calf.median=quantile(LH$calf,c(.5))
calf.ci=quantile(LH$calf,c(.025,.975))
adult.male.median=quantile(LH$adult.male,c(.5))
adult.male.ci=quantile(LH$adult.male,c(.025,.975))
lambda.healthy=list(lambda.mean=lambda.mean,lambda.median=lambda.median, lambda.ci=lambda.ci, adult.fem.median=adult.fem.median, adult.fem.ci=adult.fem.ci, yearling.fem.median=yearling.fem.median,yearling.fem.ci=yearling.fem.ci, calf.median=calf.median, calf.ci=calf.ci,adult.male.median=adult.male.median, adult.male.ci=adult.male.ci)


save(lambda.healthy, file=path("lambda_healthy.Rdata", r=root2))


LI=lambda_infected
hist(LI$lambda, breaks=60, main="B",  freq=FALSE,xlim=c(.90,1.25), family="serif", xlab = expression(paste("Population growth rate, ", lambda)))
lambda.mean=mean(LI$lambda)
lambda.median=median(LI$lambda)
lambda.sd=sd(LI$lambda)
lambda.ci=quantile(LI$lambda,c(.025,.975))
text(.96,15,paste("mean = ",signif(lambda.mean,3)),family="serif")
text(.96,12,paste("median = ",signif(lambda.median,3)),family="serif")
text(.96,9,paste("SD = ",signif(lambda.sd,3)),family="serif")
text(.96,6,paste("BCI = ",signif(lambda.ci[1],3), ", ",signif(lambda.ci[2],3)),family="serif")

adult.fem.median=quantile(LI$adult.fem,c(.5))
adult.fem.ci=quantile(LI$adult.fem,c(.025,.975))
yearling.fem.median=quantile(LI$yearling.fem,c(.5))
yearling.fem.ci=quantile(LI$yearling.fem,c(.025,.975))
calf.median=quantile(LI$calf,c(.5))
calf.ci=quantile(LI$calf,c(.025,.975))
adult.male.median=quantile(LI$adult.male,c(.5))
adult.male.ci=quantile(LI$adult.male,c(.025,.975))
lambda.infected=list(lambda.mean=lambda.mean,lambda.median=lambda.median, lambda.ci=lambda.ci, adult.fem.median=adult.fem.median, adult.fem.ci=adult.fem.ci, yearling.fem.median=yearling.fem.median,yearling.fem.ci=yearling.fem.ci, calf.median=calf.median, calf.ci=calf.ci,adult.male.median=adult.male.median, adult.male.ci=adult.male.ci)
save(lambda.infected, file=path("lambda_infected.Rdata", r=root2))




x=mean(LI$lambda.diff)
ci=quantile(LI$lambda.dif,c(.025,.975))
lambda.diff=list(x=x,ci=ci)
save(lambda.diff,file=path("lambda_diff.Rdata", r=root2))
dev.off()

demog=as.data.frame(matrix(nrow=8,ncol=5))
demog[,1]=c("Infected","","","","Healthy","","","")
demog[,2]=c("Calves","Yearling females", "Adult females", "Non-calf males","Calves","Yearling females", "Adult females", "Non-calf males")
colnames(demog)=c("Status", "Stage", "Median", ".025", ".975")
LI=lambda.infected;LH=lambda.healthy
demog[1,3:5] = c(LI$calf.median, LI$calf.ci)
demog[2,3:5] = c(LI$yearling.fem.median, LI$yearling.fem.ci)
demog[3,3:5] = c(LI$adult.fem.median, LI$adult.fem.ci)
demog[4,3:5] = c(LI$adult.male.median, LI$adult.male.ci)
demog[5,3:5] = c(LH$calf.median, LH$calf.ci)
demog[6,3:5] = c(LH$yearling.fem.median, LH$yearling.fem.ci)
demog[7,3:5] = c(LH$adult.fem.median, LH$adult.fem.ci)
demog[8,3:5] = c(LH$adult.male.median, LH$adult.male.ci)

print(xtable(demog,align=c("l","l","l","c","c","c"),digits=2,display=c("s","s","s","fg","fg","fg")), NA.string="NA", floating=FALSE,sanitize.colnames.function = function(x) {x},include.rownames=FALSE)

save(demog,file=path("Demography table.Rdata",r=root2) )

#####Total population size
####function for error bars
#Total population size=============================
pdf(file=path("Total_N.pdf",r=root2), width =6, height = 6)
y=summary(zm$N.total, quantile, c(0.0275, 0.50, 0.975))$stat

plot.BCI(x=yr,data.y=data$census$N$count.mean[range], data.bar=2*data$census$N$count.sd[range],y.low=0,y.high=6000, model.y =y[,range],ylab="Population size", xlab = "Year")
points(yr,data$removals$count$removal[range],typ="h", lwd=5, col=topo.colors(1, alpha=.3))
dev.off()




#POPULATION COMPOSITION==================
pdf(file=path("Composition.pdf",r=root2), width =7, height = 7)

par(cex.lab=1.25, mfrow=c(2,2), font.main=1, family="serif")
text.xy=c(1990,.8)

zm=One_beta.jags
yr=(1975:2010)
range=c(6:41)
#calf ratio data========
y=summary(zm$calf.ratio, quantile, c(0.0275, 0.50, 0.975))$stat

years=as.data.frame(yr);names(years)[1]="year"

y.data=merge(data$aerial$calf[,c(1,3,4)],years, by="year", all.y=TRUE)
plot_time_series(x=y.data$year,y=y[,range],y.data=y.data$mean, ybar=2*y.data$sd, y.low=0, y.high=.3,ylab="Proportion of populaion",  xlab="",main="A. Juveniles")
load(path("lambda_healthy.Rdata",r=root2))
load(path("lambda_infected.Rdata",r=root2))
abline(h=lambda.healthy$calf.ci[1], lty="dashed")
abline(h=lambda.healthy$calf.ci[2], lty="dashed")
abline(h=lambda.infected$calf.ci[1], lty="dotted")
abline(h=lambda.infected$calf.ci[2], lty="dotted")



index=range

#yearling female proportions=====
years=as.data.frame(yr);names(years)[1]="ground.year"
y.data=merge(data$ground$mu[,c(1,5)],years, by="ground.year", all.y=TRUE)
y.data=merge(y.data,data$ground$sigma[,c(1,5)], by="ground.year",all.x=TRUE); y.data[y.data$ground.year==2000 | y.data$ground.year==2001,2:3]=NA
y=summary(zm$yrl.ratio, quantile, c(0.0275, 0.50, 0.975))$stat
plot_time_series(x=y.data$ground.year,y=y[,index],y.data=y.data$yrlg_mean, ybar=2*y.data$yrlg_sd, y.low=0, y.high=.3, ylab="", xlab="", main="B. Yearling females")
abline(h=lambda.healthy$yearling.fem.ci[1], lty="dashed")
abline(h=lambda.healthy$yearling.fem.ci[2], lty="dashed")
abline(h=lambda.infected$yearling.fem.ci[1], lty="dotted")
abline(h=lambda.infected$yearling.fem.ci[2], lty="dotted")



#adult female proportions=====
years=as.data.frame(yr);names(years)[1]="ground.year"
y.data=merge(data$ground$mu[,c(1,4)],years, by="ground.year", all.y=TRUE)
y.data=merge(y.data,data$ground$sigma[,c(1,4)], by="ground.year",all.x=TRUE)
y=summary(zm$cow.ratio, quantile, c(0.0275, 0.50, 0.975))$stat
plot_time_series(x=y.data$ground.year,y=y[,index],y.data=y.data$cow_mean, ybar=2*y.data$cow_sd, y.low=0, y.high=.6, ylab="Proportion of population", xlab="Year", main="C. Adult females")
abline(h=lambda.healthy$adult.fem.ci[1], lty="dashed")
abline(h=lambda.healthy$adult.fem.ci[2], lty="dashed")
abline(h=lambda.infected$adult.fem.ci[1], lty="dotted")
abline(h=lambda.infected$adult.fem.ci[2], lty="dotted")


#yearling and adult male proportions=====
years=as.data.frame(yr);names(years)[1]="ground.year"
y.data=merge(data$ground$mu[,c(1,3)],years, by="ground.year", all.y=TRUE)
y.data=merge(y.data,data$ground$sigma[,c(1,3)], by="ground.year",all.x=TRUE)
y=summary(zm$bull.ratio, quantile, c(0.0275, 0.50, 0.975))$stat
plot_time_series(x=y.data$ground.year,y=y[,index],y.data=y.data$bull_mean, ybar=2*y.data$bull_sd, y.low=0, y.high=.6, ylab="", xlab="Year", main="D. Yearling and adult males")
abline(h=lambda.healthy$adult.male.ci[1], lty="dashed")
abline(h=lambda.healthy$adult.male.ci[2], lty="dashed")
abline(h=lambda.infected$adult.male.ci[1], lty="dotted")
abline(h=lambda.infected$adult.male.ci[2], lty="dotted")


dev.off()
	
#SEROLOGY
pdf(file=path("Serology.pdf", r=root2), height=8,width=5)
par(mfrow=c(3,1), cex.lab=1, family="serif", font.main=1)
#calves===========
a=data$sero$calf[,3]+1
b=data$sero$calf[,4]-data$sero$calf[,3]+1
mu=a/(a+b)
stdev = sqrt(a*b/((a=b^2)*(a+b+1)))
y.data=as.data.frame(cbind(data$sero$calf,mu,stdev))
#year=as.data.frame(1984:2010)
year=as.data.frame(yr)
names(year)=c("year")
y.data=merge(y.data,year,by="year",all.y=TRUE)
y=summary(zm$pv.calf, quantile, c(0.0275, 0.50, 0.975))$stat
plot_time_series(x=y.data$year,y=y[,range],y.data=y.data$mu, ybar=y.data$stdev, y.low=0, y.high=1, ylab="Proportion  sero-positive", xlab="",main="A. Juveniles")

#yearlings===========
a=data$sero$yrl[,3]+1
b=data$sero$yrl[,4]-data$sero$yrl[,3]+1
mu=a/(a+b)
stdev = sqrt(a*b/((a=b^2)*(a+b+1)))
y.data=as.data.frame(cbind(data$sero$yrl,mu,stdev))
year=as.data.frame(yr)
names(year)=c("year")
y.data=merge(y.data,year,by="year",all.y=TRUE)
y=summary(zm$pv.yrl, quantile, c(0.0275, 0.50, 0.975))$stat
plot_time_series(x=y.data$year,y=y[,range],y.data=y.data$mu, ybar=y.data$stdev, y.low=0, y.high=1, ylab="Proportion sero-positive", xlab="",main="B. Yearling females")

#cows===========
a=data$sero$cow[,3]+1
b=data$sero$cow[,4]-data$sero$cow[,3]+1
mu=a/(a+b)
stdev = sqrt(a*b/((a=b^2)*(a+b+1)))
y.data=as.data.frame(cbind(data$sero$cow,mu,stdev))
year=as.data.frame(yr)
names(year)=c("year")
y.data=merge(y.data,year,by="year",all.y=TRUE)
y=summary(zm$pv.cow, quantile, c(0.0275, 0.50, 0.975))$stat
plot_time_series(x=y.data$year,y=y[,range],y.data=y.data$mu, ybar=y.data$stdev, y.low=0, y.high=1, ylab="Proportion sero-positive", xlab="Year", main="C. Adult females")
dev.off()
serology = list(
	juv.mean = summary(zm$xpv.calf, mean)$stat,
	juv.ci = summary(zm$xpv.calf, quantile, c(.025,.975))$stat,
	yrl.mean = summary(zm$xpv.yrl, mean)$stat,
	yrl.ci = summary(zm$xpv.yrl, quantile, c(.025,.975))$stat,
	cow.mean = summary(zm$xpv.cow, mean)$stat,
	cow.ci = summary(zm$xpv.cow, quantile, c(.025,.975))$stat
)
save(serology,file=path("serology.Rdata", r=root2))



#proportion infectious
pdf(file=path("proportion infectious.pdf",r=root2), height = 7, width =4)
par(mfrow=c(2,1), family="serif")
y=summary(zm$in.per, quantile, c(0.0275, 0.50, 0.975))$stat
plot(yr,y[2,range],typ="l", xlab="", ylab="Proportion infectious", main="A", ylim=c(0,.8))
xx <-c(yr,rev(yr))
yy <- c(c(y[3,range]),rev(c(y[1,range])))
polygon(xx,yy,col=topo.colors(5,alpha=.3), border=NA)

y=summary(zm$in.per.pos, quantile, c(0.0275, 0.50, 0.975))$stat

plot(yr,y[2,range],typ="l", xlab="Year", ylab="Proportion infectious", main="B", ylim=c(0,.8))
xx <-c(yr,rev(yr))
yy <- c(c(y[3,range]),rev(c(y[1,range])))
polygon(xx,yy,col=topo.colors(5,alpha=.3), border=NA)
points(2008,.239,pch=19)
arrows(2008,.163,2008,.315, code=3, angle=90, length=.01)

#Roffe proportion infectious
a=12+1;b=26-12+1
meanp = a/(a+b)
lowp= qbeta(.025,a,b)
highp=qbeta(.975,a,b)
points(1996,meanp,pch=19)
arrows(1996,lowp,1996,highp, code=3, angle=90, length=.01)

dev.off()

x=summary(zm$in.per.pos,mean)$stat[39]
ci=summary(zm$in.per.pos,quantile,c(.025,.975))$stat[,39]
prop.infect.2008=list(x=x,ci=ci)
save(prop.infect.2008, file=path("prop_infect_2008.Rdata",r=root2))
x=summary(zm$in.per.pos,mean)$stat[27]
ci=summary(zm$in.per.pos,quantile,c(.025,.975))$stat[,27]
prop.infect.1996=list(x=x,ci=ci)
save(prop.infect.1996, file=path("prop_infect_1996.Rdata",r=root2))

prop.infect=list(
	pop.x=summary(zm$xin.per,mean)$stat,
	pop.ci=summary(zm$xin.per,quantile,c(.025,.975))$stat,
	pos.x=summary(zm$xin.per.pos,mean)$stat,
	pos.ci=summary(zm$xin.per.pos,quantile,c(.025,.975))$stat


)
save(prop.infect, file=path("prop.infect.Rdata",r=root2))

pdf(file=path("Probability of infection and Re.pdf", r=root2) ,height=8, width=5)
par(mfrow=c(2,1), family="serif")
#effective reproductive rati0
y=summary(zm$Re, quantile, c(0.0275, 0.50, 0.975))$stat
plot(yr,y[2,range],typ="l", xlab="", ylab=expression(paste("Effective reproductive ratio, ", Re[t])), main="A", ylim=c(0,4))
xx <-c(yr,rev(yr))
yy <- c(c(y[3,range]),rev(c(y[1,range])))
polygon(xx,yy,col=topo.colors(5,alpha=.3), border=NA)
abline(h=1, lty="dashed")

q_xphi=summary(zm$xphi,quantile, c(.025,.5,.975))
save(q_xphi,file="Average infection probability.Rdata")



#probability of horizontal infection aka force of infection
year=as.data.frame(yr)
names(year)="year"
y.data=as.data.frame(data$trans_data$trans)
y.data=merge(y.data,year,by="year",all.y=TRUE)
y=summary(zm$phi, quantile, c(0.0275, 0.50, 0.975))$stat
plot_time_series(x=y.data$year,y=y[,range],y.data=y.data$mean, ybar=2*y.data$sd, y.low=0, y.high=.8, ylab=expression(paste("Probability of horizontal infection, ", phi[t])), xlab="Year",main="B")

dev.off()

pdf(file=path("Removals and phi Re.pdf", r=root2) ,height=8, width=5)

plot(data$removals$count$removal[range],summary(zm$Re, quantile, c(0.50))$stat[range], xlab="", ylab=expression(paste("Probability of horizontal infection, ", phi[t])), main = "A")
plot(data$removals$count$removal[range],summary(zm$Re, quantile, c(0.50))$stat[range], xlab="Number of animals removed", ylab=expression(paste("Effective reproductive ratio, ", Re[t])), main="B")
dev.off()



##statistic for text
xphi=summary(zm$xphi,mean)$stat
ci.phi=summary(zm$xphi, quantile, c(.025,.975))
phi=list(xphi,ci.phi)
save(phi,file=path("phi_x.Rdata", r=root2 ))


#out of sample plot
pdf(file=path("Probability of infection out of sample.pdf", r=root2) ,height=6, width=6)
par(mfrow=c(1,1))
year=as.data.frame(1995:2008)
names(year)="year"
y.data=as.data.frame(data$trans_data$trans)
y.data=merge(y.data,year,by="year",all.y=TRUE)
y=summary(zm$phi, quantile, c(0.0275, 0.50, 0.975))$stat
plot_time_series(x=y.data$year,y=y[,26:39],y.data=y.data$mean, ybar=2*y.data$sd, y.low=0, y.high=1, ylab=expression(paste("Probability of Infection, ", phi)), xlab="Year",main="")
xphi=summary(zm$xphi,mean)$stat
ci.phi=summary(zm$xphi, quantile, c(.025,.975))
phi=list(xphi,ci.phi)
abline(h=.1886141, lty="dashed")
abline(h=.103, lty="dashed")
abline(h=.25797,lty="dashed")
dev.off()
save(phi,file=path("phi_x_out_of_sample.Rdata", r=root2 ))



###Probabilites of meeting objectives
###without management uncertainty
load(path("Forecasts/Without managment uncertainty/One_beta_jags.Rdata"))

z=One_beta.jags
source("Figures and Tables/Effect of management plotsII.R")

###NOTE!!! 
#the ci[]'s are calculated in the lambda plotting section, above


ecdf(z$calf.ratio.s[6,1,,])(calf.ci[1])
#calf ratios within limits

#data structure



P=array(0,dim=c(6,5,7))# 6 years, 5 managment alternatives, 7 goals
P.diff=array(0,dim=c(6,5,7))
P.ratio=array(0,dim=c(6,5,7))


for(s in 1:5){
for(t in 1:6){
		#calves in population
		#probablity of meeting objective
		upper = ecdf(z$calf.ratio.s[t,s,,])(calf.ci[2])
		lower = ecdf(z$calf.ratio.s[t,s,,])(calf.ci[1])
		P.mgmt = upper - lower
		#effect of management
		upper = ecdf(z$calf.ratio.s[t,1,,])(calf.ci[2])
		lower = ecdf(z$calf.ratio.s[t,1,,])(calf.ci[1])
		P.nothing = upper-lower
	    diff = P.mgmt - P.nothing
	    
		P[t,s,1] = P.mgmt
		P.diff[t,s,1] = diff
		ratio=P.mgmt/P.nothing
		P.ratio[t,s,1] = ratio
		
		
		#do plots for vacination example
		par(mfrow=c(2,1))
				
		plot(density(z$calf.ratio[t,2,,]),col="red", ylim=c(0,12))
		lines(density(z$calf.ratio[t,1,,]))
		abline(v=calf.ci[1])
		abline(v=calf.ci[2])

		
				
		#bulls in population
		#probablity of meeting objective
		upper = ecdf(z$bull.ratio.s[t,s,,])(adult.male.ci[2])
		lower = ecdf(z$bull.ratio.s[t,s,,])(adult.male.ci[1])
		P.mgmt = upper - lower
	
		#effect of management
		upper2 = ecdf(z$bull.ratio.s[t,1,,])(adult.male.ci[2])
		lower2 = ecdf(z$bull.ratio.s[t,1,,])(adult.male.ci[1])
		P.nothing = upper2-lower2
	    diff = P.mgmt - P.nothing
		P[t,s,2] = P.mgmt
		P.diff[t,s,2] = diff
		ratio=P.mgmt/P.nothing
		P.ratio[t,s,2] = ratio
		
		#cows in population
		#probablity of meeting objective
		upper = ecdf(z$cow.ratio.s[t,s,,])(adult.fem.ci[2])
		lower = ecdf(z$cow.ratio.s[t,s,,])(adult.fem.ci[1])
		P.mgmt = upper - lower
	
		#effect of management
		upper2 = ecdf(z$cow.ratio.s[t,1,,])(adult.fem.ci[2])
		lower2 = ecdf(z$cow.ratio.s[t,1,,])(adult.fem.ci[1])
		P.nothing = upper2-lower2
	    diff = P.mgmt - P.nothing
		P[t,s,3] = P.mgmt
		P.diff[t,s,3] = diff
		ratio=P.mgmt/P.nothing
		P.ratio[t,s,3] = ratio
		
		#population size too large
		#probablity of meeting objective <3500
		upper = 1-ecdf(z$n.s.total[t,s,,])(3500)
		P.mgmt = upper
	
		#effect of management
		upper2 = 1-ecdf(z$n.s.total[t,1,,])(3500)
		P.nothing = upper2 
	    diff = P.mgmt - P.nothing
		P[t,s,4] = P.mgmt
		P.diff[t,s,4] = diff
		ratio=P.mgmt/P.nothing
		P.ratio[t,s,4] = ratio
		
		#population size too small
		#probablity of meeting objective <3000
		upper = ecdf(z$n.s.total[t,s,,])(3000)
		P.mgmt = upper #prevent 0
	
		#effect of management
		upper2 = ecdf(z$n.s.total[t,1,,])(3000)
		P.nothing = upper2
	    diff = P.mgmt - P.nothing
		P[t,s,5] = P.mgmt
		P.diff[t,s,5] = diff
		ratio=P.mgmt/(P.nothing)
		P.ratio[t,s,5] = ratio
			

		
		
		
		#sero-prevalence
		#probablity of meeting objective
		
		lower = ecdf(z$pv.cow.s[t,s,,])(.40) + .0000000001  # prevent 0 values
		P.mgmt = lower
	
		#effect of management
		
		lower2 = ecdf(z$pv.cow.s[t,1,,])(.40) + .0000000001
		P.nothing = lower2
	    diff = P.mgmt - P.nothing
		P[t,s,6] = P.mgmt
		P.diff[t,s,6] = diff
		ratio=P.mgmt/P.nothing
		P.ratio[t,s,6] = ratio
		

		# #Change in phi
		# #probablity of meeting objective
		
		# lower = ecdf(z$delta.phi.s[t,s,,])(0)
		# P.mgmt = lower
	
		# #effect of management
		
		# lower2 = ecdf(z$delta.phi.s[t,1,,])(0)
		# P.nothing = lower2
	    # diff = P.mgmt - P.nothing
		# P[t,s,7] = P.mgmt
		# P.diff[t,s,7] = diff
		# ratio=P.mgmt/P.nothing
		# P.ratio[t,s,7] = ratio

#Cut phi in half
		#probablity of meeting objective
		
		lower = ecdf(z$phi.s[t,s,,])(summary(z$phi.s,median)$stat[1,1]/2)
		P.mgmt = lower
	
		#effect of management
		
		lower2 = ecdf(z$phi.s[t,1,,])(summary(z$phi.s,median)$stat[1,1]/2)

		P.nothing = lower2
	    diff = P.mgmt - P.nothing
		P[t,s,7] = P.mgmt
		P.diff[t,s,7] = diff
		ratio=P.mgmt/P.nothing
		P.ratio[t,s,7] = ratio


		}
	}
#If you get the error Error in sort(x) : subscript out of bounds, try stepping through each s, adding a new s each time.

	
plot(density(z$delta.phi.s[6,2,,]),xlim=c(-1,1))
lines(density(z$bull.ratio[2,1,,]), col="red")

s=as.data.frame(matrix(NA,nrow=35,ncol=8))
colnames(s)=c("Scenario", "Goal", "Year 0", "Year 1", "Year 2", "Year 3", "Year 4", "Year 5")
count=0
for(i in 1:5){  #scenario
	for(j in 1:7){ #goal
		count=count+1
		s[count,3]=signif(P[1,i,j],2)
		for(t in 2:6){ #years with treatment
		s[count,t+2] = paste(signif(P[t,i,j],2),"(", signif(P.ratio[t,i,j],2),")")

}}}
goal=c("Juveniles", "Males", "Females", "N > 3000","N < 3500","Sero-prevalence", "Pr(Transmission)")
s[,2]=rep(goal,5)
s[1,1]="No action"
s[2:7,1]= " "
s[8,1] = "Remove sero-positives"
s[9:14,1] = " "
s[15,1] = "Remove sero-negatives"
s[16:21,1] = " "
s[22,1] = "Hunting"
s[23:28,1] = " "
s[29,1] = "Vaccination"
s[30:35,1] = " "

library(xtable)

print(xtable(s[,c(1:4,6,8)]), align=c("l", "l", "l", rep("c",6)), digits=2, display=c(rep("s",3), rep("fg",6)),NA.string="NA", floating=FALSE,include.rownames=FALSE )


#print(xtable(s), align=c("l", "l", "l", rep("c",6)), digits=2, display=c(rep("s",3), rep("fg",6)),NA.string="NA", floating.environment='sidewaystable', floating=TRUE, include.rownames=FALSE )

save(s, file=path("Mgmt scenarios table.Rdata",r=root2))
#do plots for vacination example
pdf(file=path("Vacination_forecast_plots.pdf", r=root2), height=8, width=6)
par(mfrow=c(4,2))
		manage_plot(z.0=z$calf.ratio[6,1,,], z.treat=z$calf.ratio[6,5,,], lower=calf.ci[1], upper=calf.ci[2],type="between",label= "Proportion of population juvenile", panel="A - Within", xlim=c(0,.4))
		
		manage_plot(z.0=z$bull.ratio[6,1,,], z.treat=z$bull.ratio[6,5,,], lower=adult.male.ci[1], upper=adult.male.ci[2],type="between",label= "Proportion of population yearling or adult male", panel="B - Within", xlim=c(0,.80))
		
		manage_plot(z.0=z$cow.ratio[6,1,,], z.treat=z$cow.ratio[6,5,,], lower=adult.fem.ci[1], upper=adult.fem.ci[2],type="between",label= "Proportion of population adult female", panel="C - Within",xlim=c(0,.8))
		
		manage_plot(z.0=z$n.s.total[6,1,,], z.treat=z$n.s.total[6,5,,],upper=3000,type="above",label= "Population size", panel="D - Above",xlim=c(0,15000))
		manage_plot(z.0=z$n.s.total[6,1,,], z.treat=z$n.s.total[6,5,,],lower=3500,type="below",label= "Population size", panel="E - Below",xlim=c(0,15000))
		
		manage_plot(z.0=z$pv.cow.s[6,1,,], z.treat=z$pv.cow.s[6,5,,],lower=.40,type="below",label= "Seroprevalance", panel="F - Below",xlim=c(.20,1))
		
		lower = ecdf(z$phi.s[6,5,,])(summary(z$phi.s,median)$stat[1,1]/2)
		manage_plot(z.0=z$phi.s[6,1,,], z.treat=z$phi.s[6,5,,],lower=lower,type="below",label= "Probability of infection", panel="G - Below",xlim=c(0,.4))
		
dev.off()
				



# forecasts ##################
stats=array(0,dim=c(4,5,6,3))# 4 years(0,1,3,5), 5 managment alternatives, 6 states (only 1 for total, there are 2 above), 3 statistics(.025, .5, .975 quantiles)
t=c(1,2,4,6)
out.calf = summary(z$calf.ratio.s,quantile,c(.025,.5,.975))$stat 
out.bull = summary(z$bull.ratio.s,quantile,c(.025,.5,.975))$stat 
out.cow = summary(z$cow.ratio.s,quantile,c(.025,.5,.975))$stat 
out.total = summary(z$n.s.total,quantile,c(.025,.5,.975))$stat
out.prev = summary(z$pv.cow.s,quantile,c(.025,.5,.975))$stat 
out.phi =  summary(z$phi.s,quantile,c(.025,.5,.975))$stat
for(s in 1:5){
for(j in 1:4){
for(q in 1:3){
		#calves in population
		stats[j,s,1,q]=out.calf[q,t[j],s]
		stats[j,s,2,q]=out.bull[q,t[j],s]
		stats[j,s,3,q]=out.cow[q,t[j],s]
		stats[j,s,4,q]=out.total[q,t[j],s]
		stats[j,s,5,q]=out.prev[q,t[j],s]
		stats[j,s,6,q]=out.phi[q,t[j],s]
		
		
}
}
}	

#plotting  and table  for TWS talk


par(cex.lab=1.25, cex=1.25)
x=c(0:5)
plot(x, out.phi[2,,1], ylim=c(0,.25), col="red", typ="b", ylab="Probility of transmission",xlab="Years into the future")
lines(x, out.phi[2,,5], ylim=c(0,.6), col="blue", pch=19, typ="b", ylab="Probility of transmission",xlab="Years into the future")
legend(0,.05,c("No action","Vaccinate"), col =c("red","blue"), pch=c(1,19))


plot(x, out.phi[2,,1], ylim=c(0,.35), col="red", typ="b", ylab="Probility of transmission",xlab="Years into the future")
lines(x,out.phi[1,,1],col="red", typ="l", lty="dashed")
lines(x,out.phi[3,,1],col="red", typ="l", lty="dashed")
lines(x, out.phi[2,,5], ylim=c(0,.6), col="blue", pch=19, typ="b", ylab="Probility of transmission",xlab="Years into the future")
lines(x,out.phi[1,,5],col="blue", typ="l", lty="dashed")
lines(x,out.phi[3,,5],col="blue", typ="l", lty="dashed")

legend(0,.05,c("No action","Vaccinate"), col =c("red","blue"), pch=c(1,19))	
plot(density(z$phi.s[6,5,,],adjust=2),xlim=c(0,.4),col="blue", lwd=2,main="Five years in the future", xlab= "Probablity of tranmission")
lines(density(z$phi.s[6,1,,],adjust=2),col="red",lwd=2)
abline(v=.08)

p.vacine = ecdf(z$phi.s[6,5,,])(.08)
p.nothing = ecdf(z$phi.s[6,1,,])(.08)

legend(.2,10,c("No action","Vaccinate"), col =c("red","blue"), lty=c("solid", "solid"), lwd=2)

plot(density(z$delta.phi.s[6,5,,],adjust=2),xlim=c(-.4,.4),col="blue", lwd=2,main="Five years in the future", xlab= "Changen in Probablity of tranmission")
lines(density(z$delta.phi.s[6,1,,],adjust=2),col="red",lwd=2)
abline(v=0)

p.vacine.delta = ecdf(z$delta.phi.s[6,5,,])(0)
p.nothing.delta = ecdf(z$delta.phi.s[6,1,,])(0)

legend(.2,10,c("No action","Vaccinate"), col =c("red","blue"), lty=c("solid", "solid"), lwd=2)


		
		
				
		
stat.out=as.data.frame(matrix(NA,nrow=30,ncol=6))
colnames(stat.out)=c("Scenario", "Variable", "Year 0", "Year 1" ,"Year 3", "Year 5")
count=0
digits=c(2,2,2,3,2,2)
for(i in 1:5){  #scenario
	for(j in 1:6){ #state
		count=count+1
	 	for(t in 1:4){ #year
		stat.out[count,t+2] = paste(signif(stats[t,i,j,2],digits[j]),"(", signif(stats[t,i,j,1],digits[j]) , ",", signif(stats[t,i,j,3],digits[j]),")")

}}}



goal=c("Juveniles", "Males", "Adult Females", "Total Population","Sero-prevalence", "Pr(Transmission)")
stat.out[,2]=rep(goal,5)
stat.out[1,1]="No action"
stat.out[2:6,1]= " "
stat.out[7,1] = "Remove sero-positives"
stat.out[8:12,1] = " "
stat.out[13,1] = "Remove sero-negatives"
stat.out[14:18,1] = " "
stat.out[19,1] = "Hunting"
stat.out[20:24,1] = " "
stat.out[25,1] = "Vaccination"
stat.out[26:30,1] = " "

library(xtable)

colnames(stat.out)[2]="Goal"

print(xtable(stat.out), align=c("l", "l", "l", rep("c",4)), digits=2, display=c(rep("s",3), rep("fg",4)),NA.string="NA", floating=FALSE,include.rownames=FALSE )

save(stat.out, file=path("Mgmt scenarios forecast table.Rdata",r=root2))

################With management uncertainty##################
load("/Users/Tom/Documents/YNP_Bison_Model/A_final_model_runs_for_revised_manuscript/Forecasts/With management uncertainty/One_beta_jags.Rdata")
###NOTE!!! 
#the ci[]'s are calculated in the lambda plotting section, above
z=One_beta.jags

ecdf(z$calf.ratio.s[6,1,,])(calf.ci[1])
#calf ratios within limits





#data structure



P=array(0,dim=c(6,5,7))# 6 years, 5 managment alternatives, 7 goals
P.diff=array(0,dim=c(6,5,7))
P.ratio=array(0,dim=c(6,5,7))






for(s in 1:5){
for(t in 1:6){
		#calves in population
		#probablity of meeting objective
		upper = ecdf(z$calf.ratio.s[t,s,,])(calf.ci[2])
		lower = ecdf(z$calf.ratio.s[t,s,,])(calf.ci[1])
		P.mgmt = upper - lower
		#effect of management
		upper = ecdf(z$calf.ratio.s[t,1,,])(calf.ci[2])
		lower = ecdf(z$calf.ratio.s[t,1,,])(calf.ci[1])
		P.nothing = upper-lower
	    diff = P.mgmt - P.nothing

		P[t,s,1] = P.mgmt
		P.diff[t,s,1] = diff
		ratio=P.mgmt/P.nothing
		P.ratio[t,s,1] = ratio
		
		plot(density(z$calf.ratio.s))
	
		#bulls in population
		#probablity of meeting objective
		upper = ecdf(z$bull.ratio.s[t,s,,])(adult.male.ci[2])
		lower = ecdf(z$bull.ratio.s[t,s,,])(adult.male.ci[1])
		P.mgmt = upper - lower
	
		#effect of management
		upper2 = ecdf(z$bull.ratio.s[t,1,,])(adult.male.ci[2])
		lower2 = ecdf(z$bull.ratio.s[t,1,,])(adult.male.ci[1])
		P.nothing = upper2-lower2
	    diff = P.mgmt - P.nothing
		P[t,s,2] = P.mgmt
		P.diff[t,s,2] = diff
		ratio=P.mgmt/P.nothing
		P.ratio[t,s,2] = ratio
		
		#cows in population
		#probablity of meeting objective
		upper = ecdf(z$bull.ratio.s[t,s,,])(adult.fem.ci[2])
		lower = ecdf(z$bull.ratio.s[t,s,,])(adult.fem.ci[1])
		P.mgmt = upper - lower
	
		#effect of management
		upper2 = ecdf(z$bull.ratio.s[t,1,,])(adult.fem.ci[2])
		lower2 = ecdf(z$bull.ratio.s[t,1,,])(adult.fem.ci[1])
		P.nothing = upper2-lower2
	    diff = P.mgmt - P.nothing
		P[t,s,3] = P.mgmt
		P.diff[t,s,3] = diff
		ratio=P.mgmt/P.nothing
		P.ratio[t,s,3] = ratio
		
		#population size too large
		#probablity of meeting objective <3500
		upper = 1-ecdf(z$n.s.total[t,s,,])(3500)
		P.mgmt = upper
	
		#effect of management
		upper2 = 1-ecdf(z$n.s.total[t,1,,])(3500)
		P.nothing = upper2
	    diff = P.mgmt - P.nothing
		P[t,s,4] = P.mgmt
		P.diff[t,s,4] = diff
		ratio=P.mgmt/P.nothing
		P.ratio[t,s,4] = ratio
		
		#population size too small
		#probablity of meeting objective <300
		upper = ecdf(z$n.s.total[t,s,,])(3000)
		P.mgmt = upper
	
		#effect of management
		upper2 = ecdf(z$n.s.total[t,1,,])(3000)
		P.nothing = upper2
	    diff = P.mgmt - P.nothing
		P[t,s,5] = P.mgmt
		P.diff[t,s,5] = diff
		ratio=P.mgmt/P.nothing
		P.ratio[t,s,5] = ratio

		
		
		
		#sero-prevalence
		#probablity of meeting objective
		
		lower = ecdf(z$pv.cow.s[t,s,,])(.40) + .0000000001  # prevent 0 values
		P.mgmt = lower
	
		#effect of management
		
		lower2 = ecdf(z$pv.cow.s[t,1,,])(.40) + .0000000001
		P.nothing = lower2
	    diff = P.mgmt - P.nothing
		P[t,s,6] = P.mgmt
		P.diff[t,s,6] = diff
		ratio=P.mgmt/P.nothing
		P.ratio[t,s,6] = ratio
		

		# #Change in phi
		# #probablity of meeting objective
		
		# lower = ecdf(z$delta.phi.s[t,s,,])(0)
		# P.mgmt = lower
	
		# #effect of management
		
		# lower2 = ecdf(z$delta.phi.s[t,1,,])(0)
		# P.nothing = lower2
	    # diff = P.mgmt - P.nothing
		# P[t,s,7] = P.mgmt
		# P.diff[t,s,7] = diff
		# ratio=P.mgmt/P.nothing
		# P.ratio[t,s,7] = ratio

#Cut phi in half
		#probablity of meeting objective
		
		lower = ecdf(z$phi.s[t,s,,])(summary(z$phi.s,median)$stat[1,1]/2)  #(summary(z$phi.s,median)$stat[1,1]/2) is the goal determined by the value of phi in year 0.  Divide by 2 to get the 50% reduction

		P.mgmt = lower
	
		#effect of management
		
		lower2 = ecdf(z$phi.s[t,1,,])(summary(z$phi.s,median)$stat[1,1]/2)

		P.nothing = lower2
	    diff = P.mgmt - P.nothing
		P[t,s,7] = P.mgmt
		P.diff[t,s,7] = diff
		ratio=P.mgmt/P.nothing
		P.ratio[t,s,7] = ratio


		}
	}
	
	#####meet goal of reducing phi in year 5 with vaccination
	
P.mgmt =  ecdf(z$phi.s[6,5,,])(summary(z$phi.s,median)$stat[1,1]*.90)
P.nothing =  ecdf(z$phi.s[6,1,,])(summary(z$phi.s,median)$stat[1,1]*.90)
ratio=P.mgmt/P.nothing
unc2=list(P.mgmt=P.mgmt,ratio=ratio)
save(unc2, file=path("Alternative goal for phi.Rdata", r=root2))

	
	# Reduce phi
	
#If you get the error Error in sort(x) : subscript out of bounds, try stepping through each s, adding a new s each time.

	
plot(density(z$delta.phi.s[6,2,,]),xlim=c(-1,1))
lines(density(z$bull.ratio[2,1,,]), col="red")

s=as.data.frame(matrix(NA,nrow=35,ncol=8))
colnames(s)=c("Scenario", "Goal", "Year 0", "Year 1", "Year 2", "Year 3", "Year 4", "Year 5")
count=0
for(i in 1:5){  #scenario
	for(j in 1:7){ #goal
		count=count+1
		s[count,3]=signif(P[1,i,j],2)
		for(t in 2:6){ #years with treatment
		s[count,t+2] = paste(signif(P[t,i,j],2),"(", signif(P.ratio[t,i,j],2),")")

}}}
goal=c("Juveniles", "Males", "Females", "N > 3000","N < 3500","Sero-prevalence", "Pr(Transmission)")
s[,2]=rep(goal,5)
s[1,1]="No action"
s[2:7,1]= " "
s[8,1] = "Remove sero-positives"
s[9:14,1] = " "
s[15,1] = "Remove sero-negatives"
s[16:21,1] = " "
s[22,1] = "Hunting"
s[23:28,1] = " "
s[29,1] = "Vaccination"
s[30:35,1] = " "

library(xtable)



#print(xtable(s), align=c("l", "l", "l", rep("c",6)), digits=2, display=c(rep("s",3), rep("fg",6)),NA.string="NA", floating.environment='sidewaystable', floating=TRUE, include.rownames=FALSE )
s_unc=s
print(xtable(s_unc[,c(1:4,6,8)]), align=c("l", "l", "l", rep("c",6)), digits=2, display=c(rep("s",3), rep("fg",6)),NA.string="NA", floating=FALSE,include.rownames=FALSE )

save(s_unc, file=path("Mgmt scenarios table unc.Rdata",r=root2))



# forecasts ##################

stats=array(0,dim=c(4,5,6,3))# 4 years(0,1,3,5), 5 managment alternatives, 6 states (only 1 for total, there are 2 above), 3 statistics(.025, .5, .975 quantiles)
t=c(1,2,4,6)
out.calf = summary(z$calf.ratio.s,quantile,c(.025,.5,.975))$stat 
out.bull = summary(z$bull.ratio.s,quantile,c(.025,.5,.975))$stat 
out.cow = summary(z$cow.ratio.s,quantile,c(.025,.5,.975))$stat 
out.total = summary(z$n.s.total,quantile,c(.025,.5,.975))$stat
out.prev = summary(z$pv.cow.s,quantile,c(.025,.5,.975))$stat 
out.phi =  summary(z$phi.s,quantile,c(.025,.5,.975))$stat
for(s in 1:5){
for(j in 1:4){
for(q in 1:3){
		#calves in population
		stats[j,s,1,q]=out.calf[q,t[j],s]
		stats[j,s,2,q]=out.bull[q,t[j],s]
		stats[j,s,3,q]=out.cow[q,t[j],s]
		stats[j,s,4,q]=out.total[q,t[j],s]
		stats[j,s,5,q]=out.prev[q,t[j],s]
		stats[j,s,6,q]=out.phi[q,t[j],s]
		
		
}
}
}	

				
		
stat.out=as.data.frame(matrix(NA,nrow=30,ncol=6))
colnames(stat.out)=c("Scenario", "Variable", "Year 0", "Year 1" ,"Year 3", "Year 5")
count=0
digits=c(2,2,2,3,2,2)
for(i in 1:5){  #scenario
	for(j in 1:6){ #state
		count=count+1
	 	for(t in 1:4){ #year
		stat.out[count,t+2] = paste(signif(stats[t,i,j,2],digits[j]),"(", signif(stats[t,i,j,1],digits[j]) , ",", signif(stats[t,i,j,3],digits[j]),")")

}}}



goal=c("Juveniles", "Males", "Adult Females", "Total Population","Sero-prevalence", "Pr(Transmission)")
stat.out[,2]=rep(goal,5)
stat.out[1,1]="No action"
stat.out[2:6,1]= " "
stat.out[7,1] = "Remove sero-positives"
stat.out[8:12,1] = " "
stat.out[13,1] = "Remove sero-negatives"
stat.out[14:18,1] = " "
stat.out[19,1] = "Hunting"
stat.out[20:24,1] = " "
stat.out[25,1] = "Vaccination"
stat.out[26:30,1] = " "

library(xtable)
stat.out_unc=stat.out
print(xtable(stat.out_unc), align=c("l", "l", "l", rep("c",4)), digits=2, display=c(rep("s",3), rep("fg",4)),NA.string="NA", floating=FALSE,include.rownames=FALSE )

save(stat.out_unc, file=path("Mgmt scenarios forecast table unc.Rdata",r=root2))

####Table for number of animals removed, hunted, or vaccinated

r.table=matrix(nrow=4,ncol=5)
colnames(r.table)=c("Median", "2.5%", "97.5%", "Min", "Max")
rownames(r.table)=c("Remove sero-positives", "Remove sero-negatives", "Hunting", "Vaccination")
r.q=summary(z$x.r.ns, quantile,c(.025,.5,.975), na.rm=TRUE)$stat
r.min=round(summary(z$x.r.ns, min)$stat,0)
r.max=round(summary(z$x.r.ns, max)$stat,0)

r.table[,1]=r.q[2,2:5]
r.table[,2]=r.q[1,2:5]
r.table[,3]=r.q[3,2:5]
r.table[,4]=r.min[2:5]
r.table[,5]=r.max[2:5]

print(xtable(r.table, align=c("l",  rep("c",5)), digits=2, display=c(rep("s",1), rep("fg",5))),NA.string="NA", floating=FALSE,include.rownames=TRUE )
save(r.table,file=path("Removals with uncertainty.Rdata", r=root2))